{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c1b27173",
   "metadata": {},
   "source": [
    "# 机器读心术之文本挖掘与自然语言处理第6课书面作业\n",
    "学号：207402  \n",
    "\n",
    "**作业内容：**  \n",
    "1. 李航书第81页例6.1  \n",
    "假设随机变量$X$有5种可能的取值$A,B,C,D,E$，满足以下约束：  \n",
    "$$\n",
    "\\begin{align*}\n",
    "&P(A)+P(B)+P(C)+P(D)+P(E)=1 \\\\\n",
    "&P(A)+P(B)=\\frac{3}{10} \\\\\n",
    "&P(A)+P(C)=\\frac{1}{2}\n",
    "\\end{align*}\n",
    "$$\n",
    "按照最大熵的思路求$P(A),P(B),P(C),P(D),P(E)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48e3af34",
   "metadata": {},
   "source": [
    "用GIS方法求解。  \n",
    "特征函数如下：\n",
    "用0表示x=A，1表示x=B，2表示x=C，3表示x=D，4表示x=E。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a1ac6d0",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "f1&=1 \\\\\n",
    "f2&=\n",
    "\\begin{cases}\n",
    "1,\\quad x = 0\\; or\\; 1\\\\\n",
    "0, \\quad else \\\\\n",
    "\\end{cases} \\\\\n",
    "f3&=\n",
    "\\begin{cases}\n",
    "1,\\quad x= 0\\; or\\; 2\\\\\n",
    "0, \\quad else \\\\\n",
    "\\end{cases} \\\\\n",
    "f4 &=\n",
    "\\begin{cases}\n",
    "0,\\quad x= 0\\\\\n",
    "1, \\quad x=1 \\; or \\; 2 \\\\\n",
    "2, \\quad else \\\\\n",
    "\\end{cases} \\\\\n",
    "C &= \\underset{x \\in \\varepsilon}{\\operatorname{max}}\\sum_{j=1}^3 f_j(x)=3\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c7c59b2",
   "metadata": {},
   "source": [
    "初始的概率密度取$p(x)=[0.2,0.1,0.3,0.3,0.1]$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c906635c",
   "metadata": {},
   "source": [
    "实现代码如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "53f13353",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After 50 turns, training finished!\n",
      "p(x)= [0.18585902 0.11415414 0.31412154 0.19293265 0.19293265]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "class MaxEntropy:\n",
    "    def __init__(self,p_init,fp,C,delta=0.001,maxturns=10000):\n",
    "        self.p_init=p_init #表示初始随机变量的取值，第1项表示第1个取值的概率，这里即P(A)，以此类推\n",
    "        self.p_train=np.copy(p_init)\n",
    "        self.n=len(p_init) #表示随机变量有几个取值\n",
    "        self.fp=fp  #特征函数列表\n",
    "        self.m=len(fp) #有多少个特征函数\n",
    "        self.delta = delta #收敛门限\n",
    "        self.maxturns= maxturns #最大迭代次数\n",
    "        self.C= C\n",
    "    def train(self):\n",
    "        t=1\n",
    "        alpha_t=np.ones((self.m),dtype='double')\n",
    "        alpha_t1=np.ones((self.m),dtype='double')\n",
    "        #得到特征函数取值\n",
    "        f=np.zeros((self.m,self.n),dtype='double')\n",
    "        for j in range(self.m):\n",
    "            for x in range(self.n):\n",
    "                f[j][x]=self.fp[j](x)\n",
    "        p_t=np.copy(self.p_init)\n",
    "        Ef_t = np.zeros((self.m),dtype='double')\n",
    "        Ewf= np.zeros((self.m),dtype='double')\n",
    "        #计算\\tilde{E}f_j\n",
    "        for j in range(self.m):\n",
    "            for x in range(self.n):\n",
    "                Ewf[j] += self.p_init[x]*f[j][x]\n",
    "        while True:\n",
    "            #计算p^{(n)}(x)\n",
    "            for x in range(self.n):\n",
    "                p_t[x]=1.0\n",
    "                for j in range(self.m):\n",
    "                    p_t[x] *= np.power(alpha_t[j], f[j][x])\n",
    "                p_t[x] *= np.pi\n",
    "            #计算E^{(n)}f_j\n",
    "            for j in range(self.m):\n",
    "                Ef_t[j] = 0.0\n",
    "                for x in range(self.n):\n",
    "                    Ef_t[j] += p_t[x]*f[j][x]\n",
    "            #计算alpha_t1\n",
    "            for j in range(self.m):\n",
    "                alpha_t1[j]=alpha_t[j]*np.power(Ewf[j]/Ef_t[j],1./self.C)\n",
    "            t+=1\n",
    "            if np.linalg.norm(alpha_t1-alpha_t)<self.delta:\n",
    "                print(f'After {t} turns, training finished!')\n",
    "                break\n",
    "            elif t > self.maxturns:\n",
    "                print('reach max turns!')\n",
    "                break\n",
    "            else:\n",
    "                alpha_t = np.copy(alpha_t1)\n",
    "        for x in range(self.n):\n",
    "            p_t[x]=1.0\n",
    "            for j in range(self.m):\n",
    "                p_t[x] *= np.power(alpha_t[j], f[j][x])\n",
    "            p_t[x] *= np.pi\n",
    "        print('p(x)=',p_t)\n",
    "        self.p_train=np.copy(p_t)\n",
    "\n",
    "def f1(x):\n",
    "    return 1\n",
    "\n",
    "def f2(x):\n",
    "    if x == 0 or x==1:\n",
    "        return 1.0\n",
    "    else:\n",
    "        return 0.0\n",
    "\n",
    "def f3(x):\n",
    "    if x == 0 or x == 2:\n",
    "        return 1.0\n",
    "    else:\n",
    "        return 0.0\n",
    "def f4(x):\n",
    "    a=[0,1,1,2,2]\n",
    "    return a[x]\n",
    "\n",
    "p_s=np.array([0.2,0.1,0.3,0.3,0.1],dtype='double')\n",
    "fp=[f1,f2,f3,f4]\n",
    "C=3.0\n",
    "me=MaxEntropy(p_s,fp,C,0.00001)\n",
    "\n",
    "me.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adcfbec0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
